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If A is a given nonsingular matrix andA' an approximate inverse of A such that N(R) < 1, where 
R=I — AX and N is any matrix norm, then it is shown that 



1 + N(R) v ' l-N(R) 

This inequality provides the means for checking the output of a matrix inversion program. 
Methods for checking the solution of a single linear system are also discussed. 
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1. Introduction 

Although numerous excellent programs exist for the inversion of a matrix A or the solution 
of a linear system Ax = b, not enough attention has been given to the question of guaranteeing 
the accuracy of the output of such a program. In fact the problem of certifying the results of com- 
putation in general is a basic one and has been largely neglected, although it should be of pressing 
concern to any user of a high speed digital computer. 

The purpose of this article (which is mainly expository) is to point out that for matrix inversion 
there is a satisfactory solution to this question, and to derive the necessary estimates. It is the 
author's hope that this article will influence users of large machines to incorporate the procedures 
suggested (or some version of them) into their programs for matrix inversion. 

2. Matrix Norms 

The principal tool for providing such error estimates is the matrix norm. A good theoretical 
discussion of this topic may be found for example in Householder's book [3], 1 chapter 2. We 
depart somewhat from the usual definition in what follows. For our purposes a matrix norm N 
will be a real-valued function defined for all rectangular matrices over the complex numbers C 
such that for any matrices A, B 

(1) N(A) ^0 ; N(A)=0if and only if /f = 0, 

(2) N{A + B) ^N{A)+N{B), provided that A and B may be added, 
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(3) N(cA) — \c\ N(A), where c is any complex number, 

(4) N(AB) ^ N(A)N(B), provided that A and B may be multiplied in the order shown. 
Properties (1), (2), (3) are the usual properties of a distance function, but (4) is peculiar to 

matrices. 

Let A = (dij) be any mX n matrix over C. Then the most important examples of norms are 






N(A) = F(A)= \ V \aij\ 2 \ (the Euclidean or Frobenius norm), 

(6) N(A) = M(A)= vWmax \a\j\ (the maximum element norm), 

(7) N(A) = S(A)= max V |a,j| (the maximum row sum norm). 

j 

It is a simple exercise to show that these are indeed norms. The reader may wish to consult 
[4] where it is proved that (5) satisfies property (4). 

The Euclidean norm has the property that it is a unitary invariant; that is, if U, Fare unitary 
matrices of respective orders m, n, then 

F(UAV)=F(A). 
This is most easily seen from the observation that 

F(A) 2 =tr(AA*), 

where A* is the conjugate transpose of A. 

This may be used to deduce that if A is square of order n and has eigenvalues Xi, X 2 , . . ., X.», 
then 

(8) |X,| 2 +|\ 2 | 2 +. . .+\k„\^F(A)\ 

with equality if and only if A is normal (that is, if and only if A is unitarily equivalent to a diagonal 
matrix). This happens if and only if AA* = A*A. Inequality (8) is known as Schur's inequality. 

A consequence of (8) is that if X is any eigenvalue of A, then |X| ^F{A). In fact any matrix 
norm has this useful property, which we state as a theorem: 

THEOREM 1: Let A be a square matrix with arbitrary complex entries. Let X be any eigenvalue of 
A. Then if N is any matrix norm, 

(9) |A|«N(A). 

Proof: Let x be a nonzero eigenvector corresponding to X, so that Ax= Xx. Then 

\k\N(x)=N(\x) 

= N{Ax) 

*N(A)N(x), 
so that 

\M*N(A), 
since N(x) ¥= 0. This completes the proof. 
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This result provides a useful criterion for the convergence of a certain power series. Let ,4 be 
square. It is well-known that the series I + A -\- A 2 -\- . . . converges if and only if every eigenvalue 
of A is of modulus less than 1, in which case/— A is nonsingular and the sum of the series is {I — A)' 1 
(a proof may be found in [4]). Thus Theorem 1 has the useful corollary that if N(A) < 1, where N is 
any matrix norm, then I— A is nonsingular and the series /+ A + A 2 -\- . . . converges to (I—A)~ l . 
Since the convergence of this series underlies many iterative schemes used in matrix computations, 
the corollary is a valuable one. 

3. Error Bounds for Approximate Inverses 

Now suppose that A is a square nonsingular matrix, and that some computational scheme has 
been used to determine an approximate inverse^ of A. Put R = I-AX. The fact that/? is close to the 
matrix does not necessarily imply that X is close to A~ x . A simple example which illustrates this 
phenomenon is furnished by choosing 




(n 1 + n 1 — n 2 \ 



X 



where n is large. Then 



- i 
n 




R = I-AX=[ 
\0 

so that R is arbitrarily close to the matrix although X is arbitrarily far from A~ x . 

We are entitled to assume that R (which ideally should be 0) satisfies N(R) < 1 for some matrix 
norm N, since otherwise X cannot even be considered an approximation to A' 1 . We now prove the 
following, which is the main result of this article: 

THEOREM 2: Let A be a square nonsingular matrix and X an approximate inverse of A such that 
N(R) < 1 , where R = I — AX and N is some matrix norm. Then 

(10) N(XR) <N(A--X)< N(XR) ^ iXmR) 

(lV) 1 + N(R) ( ' l-N(R) l-N(R) 

PROOF: Since N(R)< 1, 1 — R is nonsinjiular and the power series I + R + R 2 + . . . converges 
to (J—R)~ l . Thus we have 

AX=I-R, 
X~ l A -'= (I-R) ' = I+R+R>+ . . . . 

Hence 

A~'-X= XR+XR1+. . . , 

N(A-'-X) « N(XR)+IV(XR*) + . . . 

« N(XR) +N(XR)N(R) + . . . 

N(XR) 
= l-N(RY 

Thus since N(XR) =£ N{X)N(R), the upper bounds given by (10) are valid. 
For the lower bound, we have 
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R = I-AX, 
XR=X-XAX 

= (A~ 1 -X)AX 

= (A-*-X)(I-R) 

= (A- 1 -X)-(A- 1 -X)R, 
that 

N(XR)=N{(A~ 1 -X)-(A-i-X)R} 
^ N(A~* -X) +N{ (A- 1 -X)R} 

^N(A-i-X)+N(A-i-X)N(R) 
= N(A~i-X) {l+N(R)}, 

N(XR) 



1 + N(R) 



\N{A-'-X). 



This completes the proof. 

It follows that for all practical purposes N(A~ l — X) may be replaced by N(XR). 

Results elaborating on Theorem 2, and applications to the evaluation of programs for matrix 
inversion, may be found in [1] and [2]. 

Theorem 2 provides the theoretical basis for inserting absolute error checks into any matrix 
inversion program, using only the information supplied to the computer and the output of the pro- 
gram. In programs which incorporate an "iterative improvement" feature there is no time loss at 
all in doing so, since the new approximantX' is computed from the odd approximantZ and the re- 
sidual R = I — AX by means of the formula 

X'=X + XR, 

so thaiXR must be calculated in any case. For matrices of moderate order, say at most 50, the time 
involved in computing XR is negligible and the sharp error bound should be used. For matrices of 
large order one multiplication may be saved by computing only R and using the weaker upper 
bound. We note parenthetically that the computation of XR must generally be carried out in double 
or higher precision, if it is to be used for iterative improvement. 

4. Error Bounds for Approximate Solutions 

The corresponding question for the solution of a linear system Ax = 6, where the output is a 
single vector which is an approximate solution of the system, does not seem capable of being 
answered so simply and completely. However, a very good relative check exists in this situation 
which is easily incorporated into any program to find x = A~ 1 b. 

Let € be the vector (1,1,..., 1) T , so that r = Ae is the vector whose coordinates are the row sums 
of A (and so available from the data). Adjoin r to the system and solve A(x, y) = (6, r), the exact 
solution of which is (x, y)= (A~ 1 b, A~ 1 r) = (A~ 1 b, e). Then the maximum deviation from 1 in the 
coordinates of the computed vector y is an excellent indication of the maximum deviation from A~ x b 
in the computed vector x. This checking procedure is both simple and effective and is readily in- 
corporated into any program which finds iL = A~ l b. 
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